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1 Introduction 



The BDDC method (Dohrmann [2003]) is one of the most advanced meth- 
i— i ods of iterative substructuring. In the case of many substructures, solving the 

coarse problem exactly becomes a bottleneck. Since the coarse problem has 
^ | the same structure as the original problem, it is straightforward to apply the 

method recursively to solve it only approximately. The two-level BDDC anal- 
ysis has been extended into three-levels in a pioneering work by Tu [2007a, b] , 
and into a general multilevel method by Mandel et al. [2008]. The methods 
for the adaptive selection of constraints for the two-level BDDC method have 

1 1 been studied in Mandel and Sousedfk [2007a], Mandel et al. [2009]. Here we 

^vq combine the two approaches into a new method preserving parallel scalability 

with increasing number of subdomains and excellent convergence properties. 

The theoretical aspects of the design of the BDDC and a closely related 
FETI-DP on irregular subdomains in the plane has been studied by Klawonn 
et al. [2008]. The authors in particular demonstrated that a proper choice 
of a certain scaling can significantly improve convergence of the methods. 
Our goal here is different. We consider only the standard stiffness scaling 
and we look for a space where the action of the BDDC preconditioner is 
defined. A combination of these two approaches, also with the proper choice of 
initial constraints (Burda et al. [2009]), would be of independent interest. The 
. £h presented algorithm has been recently extended into 3D in Sousedfk [2010]. 

^ Klawonn and Rheinbach [2007, 2008] have recently successfully developed 

5h and extensively used several inexact solvers for the FETI-DP method, and Tu 

[2008] has extended the three-level BDDC for the saddle point problems. 

All abstract spaces in this paper are finite dimensional. The dual space of 
a linear space U is denoted by U', and (-, •) is the duality pairing. 
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2 Abstract BDDC for a Model Problem 

Let C K 2 be a bounded polygonal domain, decomposed into TV nonover- 
lapping polygonal substructures J?,, i = 1,...,N, which form a conforming 
triangulation. That is, if two substructures have a nonempty intersection, 
then the intersection is a vertex, or a whole edge. Substructure vertices will 
also be called corners. Let Wj be the space of Lagrangean PI or Ql finite 
element functions with characteristic mesh size h on f2i, and which are zero 
on the boundary df2. Suppose that the nodes of the finite elements coincide 
on edges common to two substructures. Let 

ff = fix---x W N , 

and let U C W be the subspace of functions that arc continuous across the 
substructure interfaces. We wish to solve the abstract linear problem 

ueU:a(u, v) = (f,v), Vw G U, (1) 

for a given / G U' , where a is a symmetric positive semidefinitc bilinear form 
on some space W D U and positive definite on U. The form a(-, •) is called 
the energy inner product, the value of the quadratic form a (u, u) is called the 

1/2 

energy of u, and the norm ||u|| a = a (u,u) is called the energy norm. The 
operator A : U i->- JJ' associated with a is defined by 

a(u, v) — {Au, v) , Vu, v G U. 

The values of functions from W at the corners and certam averages over 
the edges will be called the coarse degrees of freedom. Let W C W be the 
space of all functions such that the values of any coarse degrees of freedom 
have a common value for all relevant substructures and vanish on df2. Define 
Ui C U C W as the subspace of all functions that are zero on all substructure 
boundaries dQi, Wa C W as the_subspace of all function such that their 
coarse degrees of freedom vanish, W n as the subspace of all functions such 
that their coarse degrees of freedom between adjacent substructures coincide 
and such that their energy is minimal. Then 

W = W A @W n , W A ± a W n . (2) 

The component of the BDDC preconditioner computed in Wn is called 
the coarse problem, cf. [Mandel et al., 2008, Algorithm 11]. Functions that 
are a-orthogonal to Ui are called discrete harmonic. In Mandel et al. [2005] 
and Mandel and Sousedfk [2007a], the analysis was done in spaces of discrete 
harmonic functions after eliminating Uj; this is not the case here, so W does 
not consist of discrete harmonic functions only. Denote by P the energy or- 
thogonal projection from U toUj. Then I — P is known as the projection onto 
the discrete harmonic functions. Finally, let £ be a projection from W onto 
U defined by taking some weighted average over the substructure interfaces. 
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Let us briefly describe the construction of the space W using the coarse 
degrees of freedom. Suppose we are given a space X and a linear operator 
C : W i-» X and define 

W = {w eW :C(I - E)w = 0} . (3) 

The values Cw will be called the local coarse degrees of freedom. To represent 
their common values, i.e. the global coarse degrees of freedom, suppose there 
is a space U c and linear operators 

Q T P :U -^U c R C :U C ^ X R:U -+W, 

such that 

CR — R c Qp. 

The operator Qp selects global coarse degrees of freedom in U c as linear 
combinations of global degrees of freedom; a global coarse degree of freedom 
is given by a row of Qp. The operator R (resp. R c ) restricts a vector of global 
(coarse) degrees of freedom into a vector of local (coarse) degrees of freedom. 
See Mandel and Sousedfk [2007a] for more details. 

2.1 Multilevel BDDC 

The substructuring components (the domains, spaces and operators) from 
the previous section will be denoted by an additional subscript i, as Q\, 
i = 1, . . . N\, etc., and called level 1. We will call the coarse problem in Wni 
the level 2 problem. It has the same finite element structure as the original 
problem (1) on level 1, so we have Ui — Wni- Level 1 substructures are level 2 
elements, level 1 coarse degrees of freedom are level 2 degrees of freedom. The 
shape functions on level 2 are the coarse basis functions in Wni, which are 
given by the conditions that the value of exactly one coarse degree of freedom 
is one and others are zero, and that they are energy minimal in W\ . Note that 
the resulting shape functions on level 2 are in general discontinuous between 
level 2 elements. Level 2 elements arc then agglomerated into nonoverlapping 
level 2 substructures, etc. Level £ elements are level i— \ substructures, and the 
level I substructures are agglomerates of level I elements. Level i substructures 
are denoted by fi\, and they are assumed to form a quasiuniform conforming 
triangulation with characteristic substructure size Hg. The degrees of freedom 
of level i elements are given by level I — 1 coarse degrees of freedom, and 
shape functions on level t are determined by^jninimization of energy on each 
level t — 1 substructure separately, so Ug = Wn,e-i- The averaging operators 
on level I, Eg : Wg -4- Ug, are defined by averaging of the values of level I 
degrees of freedom between level I substructures Q\. The space Ujg consists of 
functions in Ue that are zero on the boundaries of all level t substructures, and 
Pg : Ug — > Uig is the a— orthogonal projection in Ug onto Ujg. For convenience, 
let Q l be the original finite elements, H = h. 
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Algorithm 1 (Multilevel BDDC, Mandel et al. [2008], Algorithm 17) 

Define the preconditioner r\ G U[ i — > U\ G U\ as follows: 
for £ = 1.....L- I, 

Compute interior pre- correction on level £, 

u u eUu-.a (u ie , z u ) = (r e , z u ) , Vz H G Un- (4) 

Ge< updated residual on level £, 

rBt € C7]e, (r B£ , w £ ) = (r £ , w £ > - a (u H , ^) , V«£ G C/^. (5) 

Find the substructure correction on level £, 

wai € W 7 ^ : a (w/if, ^) = (rBi, EeZAt) , Vz^ G Wai- (6) 

Formulate the coarse problem on level £, 

w m € VK/7£ : a (wjt^, zji*) = (r B £, S^jm) , ^z m G Wjh, (7) 

If I = L — \, solve the coarse problem directly and set ui = WnL-i, 
otherwise set up the right-hand side for level £+1, 

r e+1 eW' ne , (r e+1 ,z e+1 ) = (r B£ ,E e z e+1 ) , \/z e+1 G W m = U t+1 , (8) 

end. 

for I = L - 1, . . . , 1, 

Average the approximate corrections on substructure interfaces on level I, 

UBt = E e (wAe + ue+-i) ■ (9) 
Compute the interior post- correction on level I, 

vie &U u :a (v u , z ie ) = a (u B e, z u ) , Vz /£ G U u . (10) 
Apply the combined corrections, 

ue = u u + u B e - v ie . (11) 

end. 

A condition number bound follows, cf. [Mandel et al. [2008], Lemma 20]. 
Lemma 1. If for some uj e > 1, for all £ = 1, . . . , L — 1, 

toe = sup _ Je(we), Je(we) = ll(J ~^ ) 2 " J ' lla , (12) 

w e e(I-P e )W e \\ w e\\a 

then the multilevel BDDC preconditioner satisfies k < oj — Y\k=i 
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3 Indicator of the Condition Number Bound 

As in Mandel and Sousedi'k [2007a] , we propose as an indicator of the condition 
number the maximum of the bounds from Lemma 1 computed by considering 
on each level t only one pair of adjacent substructures s and t at a time: 

w w lu = n]c~i maxwf*, tjf = sup Jf (wf) , (13) 
st t«|*e(J-p/*)M7* 

where a pair of substructures is called adjacent if they share an edge, and the 
quantities with the superscript st are defined using the domain consisting of 
the level i substructures s and t only. 

The quantity ui is called an indicator of the condition number bound. 

Let Sf be the Schur complement operator associated with the bilinear 
form a{-.-) on the space (I-Pf)Wf. The next theorem is [Mandel and 
Sousedfk, 2006, Theorem 2] written in a way suitable for our purposes. 

Theorem 2. Let a > 0, nf be the orthogonal projection onto (I — Pf)Wf , 

s f 

and I — ll t the orthogonal projection onto 

nu\\(nfsfnf + a(i-nf)). 

Then the stationary values w|\ > uf 2 > . . . and the corresponding stationary 
vectors wf k of the Rayleigh quotient Jf in (13) satisfy 

Xfwf k = u&Yfwfr (14) 
with Y e st positive definite, where 

xf = nf (i - Eff sf (i - Ef) nf, 

Yf = (nf (nfsfnf + a (i- nf)) nf + a (i- nf)) . 

The eigenvalue problem (14) is obtained by projecting the gradient of the 
Rayleigh quotient Jf (wf) onto the complement in (I — Pf )Wf of the sub- 
space, where its denominator is zero, in two steps. Both projections nf and 

— st 

n i are computed by matrix algebra, which is straightforward to implement 
numerically The projection nf projects onto null Cf (I — Ef), and I — n e 
projects onto a subspace of null Sf, which can be easily constructed com- 
putationally if a matrix Zf is given such that null Sf C range Zf . For this 
purpose, the rigid body modes are often available directly or they can be com- 
puted from the geometry of the finite element mesh. For levels I > 1, we can 
use the matrix Zf with columns consisting of coarse basis functions, because 
the span of the coarse basis functions contains the rigid body modes. In this 
way, we can reduce (14) to a symmetric eigenvalue problem, which is easier 
and more efficient to solve numerically. 
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4 Optimal Coarse Degrees of Freedom 

Writing Wf= null Cf (I - Ef) suggests how to add coarse degrees of free- 
dom to decrease the value of indicator uj. The following theorem is an anal- 
ogy to [Mandel and Sousedfk, 2006, Theorem 3]. It follows immediately from 
the standard characterization of eigenvalues as minima and maxima of the 
Rayleigh quotient on subspaces spanned by eigenvectors, applied to (14). 

Theorem 3. Suppose nf > and the coarse dof selection matrix Cf (I — Ef) 
is augmented by the rows wf k (I — Ef) T Sf (I — Ef), where wf k are the 
eigenvectors from (H). Then uf — uf nSt+1 , and ujf > ujf nSt+1 for any 

other augmentation by at most nf columns. 

In particular, if uif nSt+1 < r for all pairs of adjacent substructures s,t 

and for all levels £ = 1, . . . , L — 1, then uj < t l_1 . 

Theorem 3 allows us to guarantee that the condition number indicator 
uj < t l ~ x for a given target value r, by adding the smallest possible number 
of coarse degrees of freedom. 

The primal coarse space selection mechanism that corresponds to this 
augmentation can be explained as follows. Let us write the augmentation as 

4% = [4, k 4, k ] = <J (i Effsf (i - Ef) , 

where c| k and c\ k are blocks corresponding to substructures s and t. It should 
be noted that the matrix Ef is constructed for a pair of substructures s, t in 
such a way that, cf., e.g. [Mandel and Sousedfk, 2007b, eq. (7)], 

TjStT T^St T JT^St 

— 1 - E e , 

where Bf e and Bf are matrices known from the FETI-DP method. In par- 
ticular, the entries of Bf are +1 for substructure s and —1 for substructure t. 
This also relates our algorithm to the one from Mandel and Sousedfk [2006], 
and Mandel and Sousedfk [2007a]. Next, let us observe that, due to the appli- 
cation of J — Ef, for the two blocks of cf k it holds that c\ k ~ —c\ k , i.e., for 
the two substructures the constraint weights have the same absolute values 
and opposite sign. Hence it is sufficient to consider only one of the two blocks, 
e.g., c\ k . The augmentation of the global coarse degrees of freedom selection 
matrix [Qp.i, quj] is constructed by adding a block of k columns computed as 

nsT„sT 

Qk,e — K-i c k . 

Each column of q k ,e defines a coarse degree of freedom associated with the 
interface of level £ substructures s and t. Because i?| is a — 1 matrix, it 
means that columns in q k ,e are formed by a scattering of the entries in c s k ' . 



On Adaptive-Multilevel BDDC 7 



5 Adaptive-Multilevel BDDC in 2D 

We describe in more detail the implementation of the algorithm. It consists 
of two main steps: (i) setup, and (ii) the loop of preconditioned conjugate 
gradients with the Adaptive-Multilevel BDDC as a preconditioner. The setup 
was outlined in the previous section, and it can be summarized as follows: 

Algorithm 4 Adding of coarse degrees of freedom in order to guarantee that 
the condition number indicator Zj < t l_1 , for a given a target value t: 

for levels t = 1 : L — 1, 

Create substructures with roughly the same numbers of degrees of freedom, 
minimizing the number of "cuts" (use a graph partitioner, e.g., METIS 4-0 
(Karypis and Kumar [1998]) with weights on both, vertices and edges). 

Find a suitable set of initial constraints (corners in 2D), and set up the 
BDDC structures for the adaptive algorithm. 

for all edges En on level £, 

Compute the largest local eigenvalues and corresponding eigenvectors, until 
the first m st is found such that A^ st < t, put k = 1, . . . , m st . 

Compute the constraint weights c s k l — [c s k c\] as 

4 = w s k tT nl st (i - Ef) T sf (i - sf) n\ st , (is) 

where Il\ st is a projection constructed using the set of initial constraints. 
Take one block, e.g., c s k and keep nonzero weights for the edge Eg. 
Add to global coarse dofs selection matrix Qp,t the k columns qk,e as 

q k/ = Rfcf. (16) 

end. 

Setup the BDDC structures for level I and check size of the coarse problem: 
if small enough, call this the level L problem, factor it directly, and exit 
from the loop. 

end. 

We remark that the adaptive algorithm is significantly simpler and easier 
to implement compared to our previous algorithm from Mandel and Sousedi'k 
[2006], and Mandel and Sousedi'k [2007a]. The constraints in (15) are generated 
from the eigenvectors by the same function that evaluates the left hand side 
in (14). Then they are "torn" into two blocks, and entries of one of them, that 
correspond to a particular edge shared by substructures s and t on the level 
£, are scattered into additional columns of the matrix Qp t e- 
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The adaptive algorithm uses matrices and operators that are readily avail- 
able in our implementation of the standard BDDC method (unlike in Mandel 
et al. [2009] this time with an explicit coarse space solve) with one excep- 
tion: in order to satisfy the local partition of unity property, cf. [Mandel and 
Sousedfk, 2007b, eq. (9)], 

we need to generate locally the weight matrices Ef. 

The substructures on higher levels are then treated as (coarse) elements 
with energy minimal basis functions. However, the number of added con- 
straints is a-priori unknown. For this reason, the coarse elements must allow 
for variable number of nodes per element, and also for variable number of de- 
grees of freedom per node. It is also essential to generate a sufficient number 
of corners as initial constraints, in particular to prevent rigid body motions 
between any pair of adjacent substructures. This topic has been addressed 
several times cf., e.g., a recent contribution by Burda et al. [2009]. 

Finally, we remark that instead of performing interior pre-correction and 
post-correction on level £ = 1, cf. eqs (4)-(5) and (10)-(11), we can benefit 
from reducing the problem to interfaces in the pre-processing step. 

6 Numerical Examples and Conclusion 

The adaptive-multilevel BDDC preconditioner was implemented in Matlab 
for the 2D linear elasticity problem (with A = 1, and \i = 2) on a square 
domain discretized by finite elements with 1 182 722 degrees of freedom. The 
domain was decomposed into 2 304 subdomains on the second level and into 
9 subdomains on the third-level. Such decomposition leads to the coarsening 
ratio Hi/Hi-i — 16, with t = 1,2. In order to test the adaptive selection of 
constraints, one single edge has been jagged on both decomposition levels, see 
Fig. 1. We have computed the eigenvalues and eigenvectors of (14) by setting 
up the matrices and using standard methods for the symmetric eigenvalue 
problem in Matlab, version 7.8.0.347 (R2009a). 

In the first set of experiments, we have compared performance of the non- 
adaptive BDDC method with 2 and 3 decomposition levels. The results are 
presented in Tables 1 and 2. As expected from the theory the convergence of 
the algorithm deteriorates when additional levels are introduced. 

In the next set of experiments, we have tested the adaptive algorithm for 
the two-level BDDC. The results are summarized in Table 6. The algorithm 
performs consistently with our previous formulation in Mandel and Sousedfk 
[2007a]. The eigenvalues associated with edges between substructures clearly 
distinguish between the single problematic edge and the others (Table 3). 
Adding the coarse dofs created from the associated eigenvectors according 
to Theorem 3 decreases the value of the condition number indicator uj and 
improves convergence at the cost of increasing the number of coarse dofs. 
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Fig. 1 . The two-level decomposition into 48 x 48 (= 2304) subdomains (top) , and the 
decomposition into 9 subdomains for the three-level method (bottom); the jagged 
edge from the lower decomposition level is indicated here by a thick line. 
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Table 1. Results for non-adaptive 2-level method. Constraints are corners, or cor- 
ners and arithmetic averages over edges, denoted as c, c+e, resp. Nc is number of 
constraints, C is size of the coarse problem related to size of a subdomain problem, 
k is the condition number estimate, it is number of iterations (tol. 10~ 8 ). 

constraint Nc C n it 

c 4794 9.3 18.41 43 
c+e 13818 26.9 18.43 32 



Table 2. Results for non-adaptive 3-level method. Headings are as in Table 1. 
constraint Nc C n it 

c 4794 + 24 1.0 67.5 74 
c+e 13818 + 48 3.0 97.7 70 



Table 3. Eigenvalues of the local problems for several pairs of subdomains s and t 
on the decomposition level 1=1 (the jagged edge is between subdomains 2 and 50). 



s t A stj i 


A st ,2 


A s t,3 


A s t,4 


A s t,5 


A s t,6 




A s t,8 


1 2 


3.8 


2.4 


1.4 


1.3 


1.2 


1.1 


1.1 


1.1 


1 49 


6.0 


3.5 


2.7 


1.4 


1.3 


1.1 


1.1 


1.1 


2 3 


5.4 


2.6 


1.6 


1.3 


1.2 


1.1 


1.1 


1.1 


2 50 


24.3 


18.4 


18.3 


16.7 


16.7 


14.7 


13.5 


13.1 


3 4 


3.4 


2.4 


1.4 


1.3 


1.1 


1.1 


1.1 


1.1 


3 51 


7.4 


4.6 


3.7 


1.7 


1.4 


1.3 


1.2 


1.1 


49 50 


12.6 


5.1 


4.3 


1.9 


1.6 


1.3 


1.2 


1.2 


50 51 


8.7 


4.8 


3.9 


1.8 


1.5 


1.3 


1.2 


1.2 


50 98 


7.5 


4.6 


3.7 


1.7 


1.4 


1.3 


1.2 


1.1 



Finally, we have tested the performance of the Adaptive-Multilevel BDDC 
for the model problem with three-level decomposition (Fig. 1). Because the 
number of coarse degrees of freedom depends on an a-priori chosen value of r 
and the coarse basis functions on level i become shape basis functions on 
level £+1, the solutions of local eigenvalue problems will depend on r as well. 
This fact is illustrated by Table 4 for r = 2, and by Table 5 for r = 10 (the 
local eigenvalues for r = 3 were essentially same as for r = 2). Comparing 
the values in these two tables, we see that lower values of r result in worse 
conditioning of the local eigenvalue problems on higher decomposition level. 
This immediately gives rise to a conjecture that it might not be desirable to 
decrease the values of r arbitrarily low in order to achieve a better convergence 
of the method. On the other hand, for the model problem, comparing the 
convergence results for the two-level method (Table 6) with the three-level 
method (Table 7), we see that with the adaptive constraints we were able to 
achieve nearly the same convergence properties for the two methods. 
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Table 4. Eigenvalues of the local problems for several pairs of subdomains s, t on 
level t — 2 with r = 2 (the jagged edge is between subdomains 2 and 5). 



s t \ at ,l 


\st,2 


Ast,3 


A 3 t,4 


A 3 i,5 


A 3 t,6 


A s t,7 


Ast,8 


1 2 


16.5 


9.0 


5.4 


2.6 


2.1 


1.4 


1.3 


1.3 


1 4 


6.5 


4.7 


1.9 


1.7 


1.3 


1.2 


1.2 


1.1 


2 3 


23.1 


9.4 


4.6 


3.2 


2.1 


1.6 


1.4 


1.3 


2 5 


84.3 


61.4 


61.4 


55.9 


55.8 


49.3 


48.0 


46.9 


3 6 


13.7 


8.8 


4.4 


2.2 


1.9 


1.4 


1.3 


1.2 


4 7 


6.5 


4.7 


1.9 


1.7 


1.3 


1.2 


1.2 


1.1 


5 6 


18.9 


13.1 


11.3 


3.8 


2.6 


2.1 


1.9 


1.5 


5 8 


17.3 


12.9 


10.8 


3.6 


2.3 


2.0 


1.8 


1.4 


8 9 


13.7 


8.8 


4.4 


2.2 


1.9 


1.4 


1.3 


1.2 



Table 5. Eigenvalues of the local problems for several pairs of subdomains s, t on 
level 1 — 2 with r = 10 (the jagged edge is between subdomains 2 and 5). 



s t \ st> l 


A s i,2 


A s i,3 


A s i,4 


A s t,5 


A s t,6 


A st ,7 


A s t,8 


1 2 


7.7 


4.5 


2.7 


1.6 


1.4 


1.2 


1.2 


1.1 


1 4 


3.6 


3.0 


1.5 


1.5 


1.2 


1.2 


1.1 


1.1 


2 3 


10.9 


4.8 


2.7 


1.7 


1.5 


1.2 


1.2 


1.1 


2 5 


23.2 


17.2 


13.7 


13.7 


12.7 


12.4 


11.0 


10.9 


3 6 


6.1 


4.2 


2.5 


1.5 


1.3 


1.2 


1.1 


1.1 


4 7 


3.6 


3.0 


1.5 


1.5 


1.2 


1.2 


1.1 


1.1 


5 6 


9.8 


6.2 


4.1 


2.1 


1.6 


1.5 


1.3 


1.2 


5 8 


8.6 


5.9 


3.9 


2.0 


1.5 


1.4 


1.2 


1.2 


8 9 


6.1 


4.2 


2.5 


1.5 


1.3 


1.2 


1.1 


1.1 



Table 6. Results for the adaptive 2-level method. Headings are same as in Table 1, 
and t is the condition number target, uj is the condition number indicator. 

r Nc C U) K it 

co(=c) 4794 9.3 - 18.41 43 
10 4805 9.4 8.67 8.34 34 
3 18110 35.3 2.67 2.44 15 
2 18305 35.7 1.97 1.97 13 



Table 7. Results for the adaptive 3-level method. Headings are as in Table 6, but 
the threshold r is now used on each of the two decomposition levels and so lo < r 2 . 

r Nc C uj k it 

oo(=c) 4794 + 24 1.0 - 67.5 74 
10 4805 + 34 1.0 84.97 37.42 60 
3 18110 + 93 3.9 7.88 3.11 19 
2 18305 + 117 4.0 3.84 2.28 15 
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